Creating codebooks

The authors provide two datasets, but no further documentation - thus, we started by attempting to automatically generate some codebooks to understand the data

library(haven)
library(gt)
library(patchwork)
library(modelsummary)
## `modelsummary` 2.0.0 now uses `tinytable` as its default table-drawing
##   backend. Learn more at: https://vincentarelbundock.github.io/tinytable/
## 
## Revert to `kableExtra` for one session:
## 
##   options(modelsummary_factory_default = 'kableExtra')
##   options(modelsummary_factory_latex = 'kableExtra')
##   options(modelsummary_factory_html = 'kableExtra')
## 
## Silence this message forever:
## 
##   config_modelsummary(startup_message = FALSE)
library(lmerTest)
## Loading required package: lme4
## Loading required package: Matrix
## 
## Attaching package: 'lmerTest'
## The following object is masked from 'package:lme4':
## 
##     lmer
## The following object is masked from 'package:stats':
## 
##     step
library(broom)
library(metafor)
## Loading required package: metadat
## Loading required package: numDeriv
## 
## Loading the 'metafor' package (version 4.6-0). For an
## introduction to the package please type: help(metafor)
library(broom.mixed)
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ tidyr::expand() masks Matrix::expand()
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ✖ tidyr::pack()   masks Matrix::pack()
## ✖ tidyr::unpack() masks Matrix::unpack()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
dat1 <- read_dta("reproduction_materials/combined data OSF.dta")

Table 1 (excluding Study 6)

Likely typo in article - Table 1 states N Study 1 as 101, but data is 102, as is in note to Fig 2 (noone was excluded in that study, so changing sample size is odd) Missing condition information for 1 out of 102 - so should have been excluded

dat1 %>% count(study, sample) %>% group_by(study) %>% 
  summarise(N_total = sum(n), N_analytic = n[sample == 1]) %>% left_join(
    
    dat1 %>% group_by(study) %>% 
      summarise(
                Real_effort = any(!is.na(numbersearches)),
                Stated_effort = any(!is.na(tokens)),
                Bonus_expectations = (any(!is.na(s2_tokens_0
    )) | any(!is.na(muchlower_score))))
) %>% gt()
## Joining with `by = join_by(study)`
study N_total N_analytic Real_effort Stated_effort Bonus_expectations
1 102 102 FALSE TRUE FALSE
2 103 96 FALSE TRUE TRUE
3 68 68 TRUE FALSE TRUE
4 306 258 TRUE FALSE TRUE
5 605 518 TRUE FALSE TRUE
# Exclude excluded participants

dat1 <- dat1 %>% filter(sample == 1)

Figure 2

s2 <- dat1 %>% 
  filter(study == 2) %>% 
  select(starts_with("s2_tokens_"), cond) %>%
  pivot_longer(
    cols = starts_with("s2_tokens"), 
    names_to = "tokens", 
    values_to = "expectation",
    names_prefix = "s2_tokens_"
  ) %>% filter(tokens != "mean") %>% 
  mutate(tokens = as.numeric(tokens)) %>% 
  group_by(tokens, cond) %>% 
  summarise(mean = mean(expectation, na.rm = TRUE),
         se = sd(expectation, na.rm = TRUE) / sqrt(sum(!is.na(expectation))),
         ci_lower = mean - 1.96 * se,
         ci_upper = mean + 1.96 * se) %>% 
ggplot(aes(x = tokens, y = mean, color = haven::as_factor(cond))) +
  geom_line() + 
  scale_y_continuous(limits = c(1, 5)) +
  geom_errorbar(aes(ymin = ci_lower, ymax = ci_upper), width = 0.2) +
  labs(x = "Tokens contributed", y = "Subjective likelihood of bonus", col = "") +
  scale_color_manual(values = c("Control" = "blue", "Advantage" = "green", "Disadvantage" = "orange")) +
  theme_minimal()
## `summarise()` has grouped output by 'tokens'. You can override using the
## `.groups` argument.
s3_5 <- dat1 %>% 
  filter(study > 2) %>% 
  select(matches("_score"), cond) %>%
  pivot_longer(
    cols = matches("_score"), 
    names_to = c("performance", NA), 
    values_to = "expectation",
        names_sep = "_"
  ) %>% 
  mutate(performance = fct_recode(performance,
                                  "Much Lower" = "muchlower",
                                  "Lower" = "lower",
                                  "About Same" = "same", 
                                  "Higher" = "higher",
                                  "Much Higher" = "muchhigher") %>% 
           factor(levels = c("Much Lower", "Lower", "About Same", "Higher", "Much Higher"))) %>%
  group_by(performance, cond) %>% 
  summarise(mean = mean(expectation, na.rm = TRUE),
         se = sd(expectation, na.rm = TRUE) / sqrt(sum(!is.na(expectation))),
         ci_lower = mean - 1.96 * se,
         ci_upper = mean + 1.96 * se) %>% 
ggplot(aes(x = performance, y = mean, color = haven::as_factor(cond), group = haven::as_factor(cond), )) +
  geom_line() + 
    scale_y_continuous(limits = c(0, 100)) +
  geom_errorbar(aes(ymin = ci_lower, ymax = ci_upper), width = 0.2) +
  labs(x = "Performance relative to other employee(s)", y = "Subjective likelihood of bonus (%)", col = "") +
  scale_color_manual(values = c("Control" = "blue", "Advantage" = "green", "Disadvantage" = "orange")) +
  theme_minimal()
## `summarise()` has grouped output by 'performance'. You can override using the
## `.groups` argument.
p1 <- s2 / s3_5

p1

Control group mean misreported in text, is mean of disadvantaged group

Means

s2_means <- dat1 %>% 
  filter(study == 2) %>%
  rowid_to_column() %>% 
  select(-tokens) %>% 
  pivot_longer(
    cols = starts_with("s2_tokens"), 
    names_to = "tokens", 
    values_to = "expectation",
    names_prefix = "s2_tokens_"
  ) %>% 
  group_by(rowid) %>% 
  summarise(cond = as_factor(cond)[1],
            m_effort = mean(as.numeric(tokens), na.rm = TRUE), # Should be fixed, just code check
            m_expectation = mean(expectation, na.rm = TRUE),
            sd_expectation = sd(expectation, na.rm = TRUE)
) 
## Warning: There were 96 warnings in `summarise()`.
## The first warning was:
## ℹ In argument: `m_effort = mean(as.numeric(tokens), na.rm = TRUE)`.
## ℹ In group 1: `rowid = 1`.
## Caused by warning in `mean()`:
## ! NAs introduced by coercion
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 95 remaining warnings.
mod_s2_avg <- lm(m_expectation ~ cond, s2_means)

modelsummary(mod_s2_avg, estimate = "{estimate} [{conf.low}, {conf.high}]",
             statistic = "({std.error}), p = {p.value}",
             gof_map = NA,
             notes = as.character(glue::glue("df = {modelsummary::glance(mod_s2_avg)$df.residual}")),
             title = "Average expectation in Study 2")
tinytable_l28mme817efhw05euxa1
Average expectation in Study 2
(1)
df = 93
(Intercept) 2.961 [2.776, 3.146]
(0.093), p = <0.001
condAdvantage 0.426 [0.172, 0.680]
(0.128), p = 0.001
condDisadvantage -0.586 [-0.844, -0.328]
(0.130), p = <0.001
s3_5_means <- dat1 %>% 
  filter(study > 2) %>%
  rowid_to_column() %>% 
  select(-tokens) %>% 
  pivot_longer(
    cols = matches("_score"), 
    names_to = c("performance", NA), 
    values_to = "expectation",
        names_sep = "_"
  ) %>% 
  group_by(rowid) %>% 
  summarise(study = study[1],
            cond = as_factor(cond)[1],
            effort_ratings = n(), # Should be fixed, just code check
            m_expectation = mean(expectation, na.rm = TRUE),
            sd_expectation = sd(expectation, na.rm = TRUE)
) 

mod_s3_5_avg <- lmer(m_expectation ~ cond  + (1 | study), s3_5_means)

modelsummary(mod_s3_5_avg,
             estimate = "{estimate} [{conf.low}, {conf.high}]",
             statistic = "({std.error}), p = {p.value}",
             title = "Average expectation in pooled studies 3-5")
tinytable_7rp384wkw4uuowymdy01
Average expectation in pooled studies 3-5
(1)
(Intercept) 49.363 [44.396, 54.330]
(2.531), p = <0.001
condAdvantage 15.203 [13.007, 17.399]
(1.119), p = <0.001
condDisadvantage -14.225 [-16.413, -12.037]
(1.115), p = <0.001
SD (Intercept study) 4.107
, p =
SD (Observations) 13.301
, p =
Num.Obs. 844
R2 Marg. 0.423
R2 Cond. 0.473
AIC 6772.1
BIC 6795.8
ICC 0.1
RMSE 13.26

Slopes

s2_expectations <- dat1 %>% 
  filter(study == 2) %>%
  rowid_to_column() %>% 
  select(-tokens) %>% 
  pivot_longer(
    cols = starts_with("s2_tokens"), 
    names_to = "tokens", 
    values_to = "expectation",
    names_prefix = "s2_tokens_"
  ) %>% filter(tokens != "mean") %>% 
  mutate(tokens = as.numeric(tokens), cond = as_factor(cond))

s2_expectations <- s2_expectations %>% group_by(cond, rowid) %>% 
  mutate(delta = expectation - lag(expectation)) %>% 
  summarise(delta = mean(delta, na.rm = TRUE))
## `summarise()` has grouped output by 'cond'. You can override using the
## `.groups` argument.
mod_s2_slope <- lm(delta ~ cond, s2_expectations) %>% summary()

modelsummary(mod_s2_slope, estimate = "{estimate} [{conf.low}, {conf.high}]",
             statistic = "({std.error}), p = {p.value}",
             gof_map = NA,
             notes = as.character(glue::glue("df = {modelsummary::glance(mod_s2_slope)$df.residual}")),
             title = "Expected reward sensitivity in Study 2")
tinytable_h1yimrxa5jam6ka81r5o
Expected reward sensitivity in Study 2
(1)
df = 93
(Intercept) 0.753 [0.660, 0.847]
(0.047), p = <0.001
condAdvantage -0.165 [-0.294, -0.037]
(0.065), p = 0.012
condDisadvantage -0.253 [-0.384, -0.123]
(0.066), p = <0.001
s3_5_expectations <- dat1 %>% 
  filter(study > 2) %>%
  rowid_to_column() %>% 
  select(-tokens) %>% 
  pivot_longer(
    cols = matches("_score"), 
    names_to = c("performance", NA), 
    values_to = "expectation",
        names_sep = "_"
  ) %>% mutate(performance = fct_recode(performance,
                                  "Much Lower" = "muchlower",
                                  "Lower" = "lower",
                                  "About Same" = "same", 
                                  "Higher" = "higher",
                                  "Much Higher" = "muchhigher") %>% 
           factor(levels = c("Much Lower", "Lower", "About Same", "Higher", "Much Higher")),
           cond = as_factor(cond)) %>% 
  group_by(study, cond, rowid) %>% 
    mutate(delta = expectation - lag(expectation)) %>% 
  summarise(delta = mean(delta, na.rm = TRUE))
## `summarise()` has grouped output by 'study', 'cond'. You can override using the
## `.groups` argument.
mod_s3_5_delta <- lmer(delta ~ cond  + (1 | study), s3_5_expectations)
## boundary (singular) fit: see help('isSingular')
modelsummary(mod_s3_5_delta,
             estimate = "{estimate} [{conf.low}, {conf.high}]",
             statistic = "({std.error}), p = {p.value}",
             title = "Expected reward sensitivity in pooled studies 3-5")
tinytable_sco7nb3uh7h213j81zk7
Expected reward sensitivity in pooled studies 3-5
(1)
(Intercept) 22.178 [21.426, 22.930]
(0.383), p = <0.001
condAdvantage -6.124 [-7.205, -5.044]
(0.551), p = <0.001
condDisadvantage -7.995 [-9.073, -6.918]
(0.549), p = <0.001
SD (Intercept study) 0.000
, p =
SD (Observations) 6.545
, p =
Num.Obs. 843
R2 Marg. 0.217
AIC 5567.0
BIC 5590.7
RMSE 6.53

Figure 2 (& table)

# Calculate standardised effects
summaries <- dat1 %>% 
  filter(!is.na(cond)) %>% 
  group_by(study) %>% 
  mutate(z_effort = (coalesce(log(numbersearches), tokens) - mean(coalesce(log(numbersearches), tokens), na.rm = TRUE))/sd(coalesce(log(numbersearches), tokens), na.rm = TRUE)) %>% 
  group_by(study, cond) %>% 
  summarise(cond = as_factor(cond)[1],
            m_effort = mean(z_effort, na.rm = TRUE),
            sd_effort = sd(z_effort, na.rm = TRUE),
            N = n(),
            se = sd_effort / sqrt(N), 
            .groups = "drop") %>% 
  mutate(conf.low = m_effort - 1.96 * se,
         conf.high = m_effort + 1.96 * se,
         vi = se^2)

# Add p-values per study
p_values <- dat1 %>%
  filter(!is.na(cond)) %>%
  group_by(study) %>%
  mutate(z_effort = (coalesce(log(numbersearches), tokens) - 
                      mean(coalesce(log(numbersearches), tokens), na.rm = TRUE)) / 
                      sd(coalesce(log(numbersearches), tokens), na.rm = TRUE),
         cond = cond %>% as_factor() %>% relevel(ref = "Control")) %>%
  group_by(study) %>%
  do({
        data <- .
    model <- lm(z_effort ~ cond, data = data)
    tidy(model) %>%
      filter(term != "(Intercept)") %>% # Remove the intercept term
      mutate(study = unique(data$study))
  }) %>%
  transmute(study, cond = term %>% str_remove("cond"), p_value = p.value) %>%
  ungroup()

summaries <- summaries %>% left_join(p_values)
## Joining with `by = join_by(study, cond)`
# Run meta-analyses

## Model without an intercept to get estimates and standard errors for each condition
model_no_intercept <- rma.mv(yi = m_effort, V = vi, mods = ~ factor(cond) - 1, random = ~ 1 | study, data = summaries)
estimates_no_intercept <- coef(summary(model_no_intercept))[, "estimate"]
se_no_intercept <- coef(summary(model_no_intercept))[, "se"]

## Model with condition as a moderator to get p-values comparing each condition to control
summaries$cond <- relevel(factor(summaries$cond), ref = "Control")
model_with_intercept <- rma.mv(yi = m_effort, V = vi, mods = ~ cond, random = ~ 1 | study, data = summaries)
summary_with_intercept <- summary(model_with_intercept)

## Extract p-values from the model with intercept
p_values <- coef(summary_with_intercept)[-1, "pval"]
p_values <- c(NA, p_values) # Control comparison is NA

meta_results <- tibble(
  study = "meta",
  cond = levels(factor(summaries$cond)),
  m_effort = estimates_no_intercept,
  se = se_no_intercept,
  p_value = p_values,
  conf.low = m_effort - 1.96 * se,
  conf.high = m_effort + 1.96 * se
)

fig2_res <- meta_results %>% bind_rows(summaries %>% mutate(study = as.character(study)), .)

# Function to calculate t-test p-value using summary statistics
fig2_res %>% select(-vi) %>% 
gt() %>% 
  fmt_number(columns = c(m_effort, conf.low, conf.high, se, sd_effort), decimals = 3) %>% 
  fmt(
    columns = p_value,
    fns = \(p) timesaveR::fmt_p(p, equal_sign = FALSE)
  )
study cond m_effort sd_effort N se conf.low conf.high p_value
1 Control 0.437 0.971 34 0.167 0.111 0.764 NA
1 Advantage 0.026 0.982 32 0.174 −0.314 0.366 .079
1 Disadvantage −0.448 0.866 35 0.146 −0.735 −0.162 < .001
2 Control 0.365 0.873 30 0.159 0.053 0.678 NA
2 Advantage 0.092 0.933 34 0.160 −0.222 0.406 .256
2 Disadvantage −0.440 1.042 32 0.184 −0.801 −0.079 .001
3 Control 0.174 1.073 23 0.224 −0.265 0.612 NA
3 Advantage 0.013 1.008 22 0.215 −0.408 0.434 .593
3 Disadvantage −0.186 0.925 23 0.193 −0.564 0.192 .228
4 Control 0.271 0.842 91 0.088 0.098 0.444 NA
4 Advantage −0.025 0.942 85 0.102 −0.226 0.175 .046
4 Disadvantage −0.274 1.143 82 0.126 −0.522 −0.027 < .001
5 Control 0.152 0.895 178 0.067 0.020 0.283 NA
5 Advantage 0.029 0.980 167 0.076 −0.120 0.177 .250
5 Disadvantage −0.184 1.093 173 0.083 −0.347 −0.021 .002
meta Control 0.020 NA NA 0.052 −0.083 0.123 NA
meta Advantage 0.229 NA NA 0.047 0.136 0.322 .003
meta Disadvantage −0.266 NA NA 0.057 −0.378 −0.155 < .001
# Convert p-values to a readable format for annotations
fig2_res <- fig2_res %>%
  mutate(p_label = ifelse(is.na(p_value), "", paste0("p ", timesaveR::fmt_p(p_value, equal_sign = TRUE, digits = 3))))

# Create the plot using ggplot2
# Reorder the conditions
fig2_res$cond <- factor(fig2_res$cond, levels = c("Control", "Advantage", "Disadvantage"))

# Create the plot using ggplot2
ggplot(fig2_res, aes(x = as.factor(study), y = m_effort, color = cond, group  = cond)) +
  geom_point(position = position_dodge(width = 0.7), size = 3, shape = 16) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.2, position = position_dodge(width = 0.7)) +
  geom_text(aes(label = p_label, y = conf.low), position = position_dodge(width = 0.7), vjust = 1.4, size = 3, color = "black") +
  labs(x = "Study", y = "Standardized mean", color = NULL) +
  theme_minimal() +
  theme(legend.position = "bottom", legend.title = element_blank(), axis.text.x = element_text(angle = 45, hjust = 1)) +
  scale_color_manual(values = c("Control" = "blue", "Advantage" = "darkgreen", "Disadvantage" = "darkorange")) +
  scale_x_discrete(labels = c("1" = "Study 1", "2" = "Study 2", "3" = "Study 3", "4" = "Study 4",
                              "5" = "Study 5", "meta" = "Mini meta-analysis")) +
  guides(color = guide_legend(override.aes = list(shape = 16))) +
  expand_limits(y = c(min(fig2_res$conf.low) - 0.1, max(fig2_res$conf.high)))

# Mediation

# Study 2
med_data <- dat1 %>% 
  filter(!is.na(cond)) %>% 
  group_by(study) %>% 
  mutate(z_perf_slope = coalesce(scale(s2_performance_slope) %>% as.numeric(),
                                 scale(performance_slope) %>% as.numeric()),
         z_numbersearches = scale(numbersearches) %>% as.numeric(),
         z_tokens = scale(tokens) %>% as.numeric(),
         cond = as_factor(cond)) %>% 
  ungroup() %>% 
  filter(!is.na(z_perf_slope)) %>% 
  mutate(timesaveR::dummy_code(cond))

Regression results for Study 2

slope_mod_s2 <- med_data %>% filter(study == 2) %>% lm(z_tokens ~ z_perf_slope, .)

slope_mod_s2 %>% tidy(conf.int = TRUE) %>% 
  mutate(df = slope_mod_s2$df.residual) %>% 
  filter(term != "(Intercept)") %>% 
  glimpse()
## Rows: 1
## Columns: 8
## $ term      <chr> "z_perf_slope"
## $ estimate  <dbl> 0.3884817
## $ std.error <dbl> 0.09504098
## $ statistic <dbl> 4.087518
## $ p.value   <dbl> 9.181902e-05
## $ conf.low  <dbl> 0.1997756
## $ conf.high <dbl> 0.5771878
## $ df        <int> 94

Regression results for Study 3-5

library(lmerTest)
slope_mod_s3_5 <- med_data %>% filter(study > 2) %>% 
  lmer(z_numbersearches ~ z_perf_slope + (1 | study), .)
## boundary (singular) fit: see help('isSingular')
slope_mod_s3_5 %>% summary()
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: z_numbersearches ~ z_perf_slope + (1 | study)
##    Data: .
## 
## REML criterion at convergence: 2375.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.6412 -0.7776 -0.2550  0.6303  2.7683 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  study    (Intercept) 0.0000   0.0000  
##  Residual             0.9714   0.9856  
## Number of obs: 843, groups:  study, 3
## 
## Fixed effects:
##               Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)  1.167e-03  3.395e-02 8.410e+02   0.034    0.973    
## z_perf_slope 1.659e-01  3.401e-02 8.410e+02   4.877 1.29e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## z_perf_slop 0.000 
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')

Mediation results

Note that N in Fig 3 is reported as 844, so that it corresponds to these results - even though it is not stated that the Fig only contains 3-5 rather than 2-5. Pattern for S2 looks very similar, though effects are larger (in line with regression results).

library(lavaan)
## This is lavaan 0.6-17
## lavaan is FREE software! Please report any bugs.
med_data <- med_data %>% mutate(z_effort = coalesce(z_tokens, z_numbersearches))

mod_full <- ("
  # Path from condition (cond) to mediator (z_perf_slope)
  z_perf_slope ~ a1 * cond_advantage + a2 * cond_disadvantage

  # Path from condition (cond) to outcome (z_effort)
  z_effort ~ Advantage * cond_advantage + Disadvantage * cond_disadvantage + b * z_perf_slope
")

mod_reduced <- ("
  z_effort ~ Advantage * cond_advantage + Disadvantage * cond_disadvantage
")

med_res <- sem(mod_full, data = med_data %>% filter(study > 2), cluster = "study") %>%
  tidy(conf.int = TRUE) %>%
  filter(str_detect(label, "Advantage|Disadvantage")) %>%
  transmute(model = "full", cond = label, study = "3-5", estimate, conf.low, conf.high, p.value) %>%
  bind_rows(
    sem(mod_reduced, data = med_data %>% filter(study > 2), cluster = "study") %>% tidy(conf.int = TRUE) %>% filter(str_detect(label, "Advantage|Disadvantage")) %>% transmute(model = "reduced", cond = label, study = "3-5", estimate, conf.low, conf.high, p.value)
  ) %>%
  bind_rows(
    sem(mod_full, data = med_data %>% filter(study == 2)) %>% tidy(conf.int = TRUE) %>% filter(str_detect(label, "Advantage|Disadvantage")) %>% transmute(model = "full", cond = label, study = "2", estimate, conf.low, conf.high, p.value) %>% bind_rows(
      sem(mod_reduced, data = med_data %>% filter(study == 2)) %>% tidy(conf.int = TRUE) %>% filter(str_detect(label, "Advantage|Disadvantage")) %>% transmute(model = "reduced", cond = label, study = "2", estimate, conf.low, conf.high, p.value)
    )
  ) %>%
  mutate(p_label = paste0("p ", timesaveR::fmt_p(p.value, equal_sign = TRUE, digits = 3)))
## Warning in lav_model_vcov(lavmodel = lavmodel, lavsamplestats = lavsamplestats, : lavaan WARNING:
##     The variance-covariance matrix of the estimated parameters (vcov)
##     does not appear to be positive definite! The smallest eigenvalue
##     (= -2.036678e-18) is smaller than zero. This may be a symptom that
##     the model is not identified.
## Warning in lav_model_vcov(lavmodel = lavmodel, lavsamplestats = lavsamplestats, : lavaan WARNING:
##     The variance-covariance matrix of the estimated parameters (vcov)
##     does not appear to be positive definite! The smallest eigenvalue
##     (= 6.488542e-20) is close to zero. This may be a symptom that the
##     model is not identified.
med_res %>%
  gt() %>%
  fmt_number(columns = c(estimate, conf.low, conf.high), decimals = 3) %>%
  fmt_number(columns = p.value, decimals = 4)
model cond study estimate conf.low conf.high p.value p_label
full Advantage 3-5 −0.035 −0.099 0.029 0.2859 p = .286
full Disadvantage 3-5 −0.154 −0.215 −0.092 0.0000 p < .001
reduced Advantage 3-5 −0.151 −0.306 0.004 0.0563 p = .056
reduced Disadvantage 3-5 −0.306 −0.458 −0.153 0.0001 p < .001
full Advantage 2 −0.083 −0.536 0.369 0.7181 p = .718
full Disadvantage 2 −0.514 −0.993 −0.036 0.0351 p = .035
reduced Advantage 2 −0.273 −0.734 0.188 0.2455 p = .245
reduced Disadvantage 2 −0.805 −1.273 −0.338 0.0007 p < .001
med_res %>%
  group_by(study, cond) %>%
  summarise(
    share_mediated = 1 - (estimate[1] / estimate[2]), order = paste(model[1], model[2]),
    .groups = "drop_last"
  ) %>%
  dplyr::select(-order) %>%
  gt() %>%
  fmt_percent(columns = share_mediated, decimals = 1)
cond share_mediated
2
Advantage 69.4%
Disadvantage 36.1%
3-5
Advantage 76.8%
Disadvantage 49.7%
med_res %>%
  mutate(model = factor(model, levels = c("reduced", "full"))) %>%
  ggplot(aes(x = cond, y = estimate, group = model, shape = model, color = model)) +
  geom_point(size = 3, position = position_dodge(width = 0.4)) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high),
    width = 0.2, position = position_dodge(width = 0.4)
  ) +
  geom_text(aes(y = conf.low - 0.05, label = p_label),
    color = "black",
    position = position_dodge(width = 0.4),
    vjust = 1, size = 3
  ) +
  facet_wrap(~study, ncol = 2, labeller = labeller(study = c("2" = "Study 2", "3-5" = "Study 3-5"))) +
  scale_shape_manual(values = c(16, 1), labels = c("Reduced model", "Full model")) +
  scale_color_manual(values = c("black", "black"), labels = c("Reduced model", "Full model")) +
  labs(x = "", y = "Standardized effect", shape = "", color = "") +
  scale_y_continuous(limits = c(-1.4, 0.5), breaks = round(seq(-1.4, 0.4, by = 0.2), 2)) +
  theme_minimal() +
  theme(
    legend.position = "bottom",
    axis.title.y = element_text(margin = margin(r = 10)),
    panel.grid.major.x = element_blank(),
    panel.grid.minor = element_blank()
  ) +
  geom_hline(yintercept = 0, linetype = "solid", color = "black")